#Replication Script - Stein and Cantin 2021

#Analysis run using R version 4.0.5

#Necessary Packages

library(tidyverse)
library(dplyr)
library(MASS)
library(readxl)
library(margins)
library(xlsx)
library(jtools)
library(sjPlot)

#Data Import

data <- Stein_Cantin_2021
data2 <- data %>% filter(Non_Lethal_Sponsorship == 0)
data3 <- data %>% filter(Lethal_Sponsorship == 0)
data4 <- Stein_Cantin_2021
data4$State_Sponsorship <- as.factor(data4$State_Sponsorship)
data5 <- data %>% filter(Location != "Syria")

#Main Models - Included in the Main Text

Model1 <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

Model2 <- glm.nb(Inter_Rebel_Fighting ~ Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data2, maxit = 800)

Model3 <- glm.nb(Inter_Rebel_Fighting ~ Non_Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data3, maxit = 800)

Model4 <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

#Plots - Included in the Main Text

Model1Plot <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data4, maxit = 800)
Model4Plot <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data4, maxit = 800)

plot_model(Model1Plot, type = "pred", terms = c("State_Sponsorship"),
           cat.geom = "line",
           jitter = .2,
           rug = TRUE,
           int.width = 0.95,
           colors = c("grey20","grey50","grey80"),
           line.thickness = 0.6,
           point.size = 0.5,
           point.alpha = 0.2,
           point.color = "black",
           cat.pred.point.size = 1.5)

plot_model(Model4Plot, type = "pred", terms = c("State_Sponsorship", "Religious"), mdrt.values = "all",
           cat.geom = "line",
           jitter = .2,
           rug = TRUE,
           int.width = 0.95,
           colors = c("grey20","grey50","grey80"),
           line.thickness = 0.6,
           point.size = 0.5,
           point.alpha = 0.2,
           point.color = "black",
           cat.pred.point.size = 1.5)

#Robustness Checks - Included in the Online Appendix

Model1a <- glm.nb(Inter_Rebel_Fighting_2 ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

Model4a <- glm.nb(Inter_Rebel_Fighting_2 ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

Model1b <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln + Conflict_Intensity_ln, data=data, maxit = 800)

Model4b <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln + Conflict_Intensity_ln, data=data, maxit = 800)

Model1c <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship_Alleged + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

Model4c <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship_Alleged*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)

Model1d <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data5, maxit = 800)

Model4d <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data5, maxit = 800)

#Model Comparison - Included in the Online Appendix

ModelA <- glm.nb(Inter_Rebel_Fighting ~ Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)
Model1 <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)
Model4 <- glm.nb(Inter_Rebel_Fighting ~ State_Sponsorship*Religious + Command + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data, maxit = 800)
anova(ModelA, Model1, Model4, test = "Chisq")

ModelB <- glm.nb(Inter_Rebel_Fighting ~ Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data2, maxit = 800)
Model2 <- glm.nb(Inter_Rebel_Fighting ~ Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data2, maxit = 800)
anova(ModelB, Model2, test = "Chisq")

ModelC <- glm.nb(Inter_Rebel_Fighting ~ Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data2, maxit = 800)
Model3 <- glm.nb(Inter_Rebel_Fighting ~ Non_Lethal_Sponsorship + Command + Religious + Ethnic + Rebel_Number + Endowments + GDP_Per_Capita_ln, data=data3, maxit = 800)
anova(ModelC, Model3, test = "Chisq")


